%{

EchochG analyzer written as sub-function for GUI app
Inputs come in as structure, expecting as examples:
in.W = W;       % Time waveform
in.fs = 3e4;    % sampling frequency
in.lb = 74;     % Start index for max/min search
in.rb = 241;    % Stop index for max/min search
in.fp1 = 3000;  % Lowpass frequency for analysis
in.fp2 = 470;   % Highpass/Lowpass frequency for analysis
in.mp = 0.02; % minimum prominence

Returns a bunch of stuff in o:
[o] = GUIsub_ECochG(in)

Rich Herr.  Mass Eye and Ear Institute
Jun2 24

%}


function [o,o2] = GUIsub_ECochG(in)

%filtering
    Wbp = highpass(lowpass(in.W,in.fp1,in.fs,'ImpulseResponse','iir','Steepness',.95),...
                             in.fp2,in.fs,'ImpulseResponse','iir','Steepness',.95);
    Wlp_3K  =  (lowpass(in.W,in.fp1,in.fs,'ImpulseResponse','iir','Steepness',.95));
    Wlp_470 = (lowpass(in.W,in.fp2,in.fs,'ImpulseResponse','iir','Steepness',.95));

    %t = (0:length(in.W)-1)*1000/in.fs-2; % MAKE THIS FLEXIBLE BASED ON INPUT TIME VECTOR, NOT HARD-CODED!
%% max: BP [0.47 3 kHz] waveform
    Wbp_max = islocalmax(Wbp(in.lb:in.rb),'MinProminence',in.mp,'MaxNumExtrema',4, 'MinSeparation',15);
    Lind_max = find(Wbp_max)+in.lb-1;
    lat_max = in.Wtime(Lind_max);
    amp_max_bp = Wbp(Lind_max);

%% max: LP < 0.47 kHz waveform
    Wmax_470 = islocalmax(Wlp_470(in.lb:in.rb),'MaxNumExtrema',3);
    Lind_max_470 = find(Wmax_470, 1, 'first')+in.lb-1;
    Lind_max_470_last = find(Wmax_470, 1, 'last')+in.lb-1;
    lat_peak_470 = in.Wtime(Lind_max_470);
    lat_peak_470_last = in.Wtime(Lind_max_470_last);
    amp_peak_470 = Wlp_470(Lind_max_470);
    amp_peak_470_last = Wlp_470(Lind_max_470_last);
    
%% min: SP and BP [0.47 3 kHz] waveform
    Wmin1 = islocalmin(Wbp(in.lb-15:Lind_max(1)),'MinProminence',in.mp-0.015,'MaxNumExtrema',3);
    Lind_min1 = find(Wmin1,1, 'last')+in.lb-15-1;
    Wmin2 = islocalmin(Wbp(Lind_max(1):in.rb-24),'MinProminence',in.mp-0.01,'MaxNumExtrema',3, 'MinSeparation', 15 );
    Lind_min2 = find(Wmin2)+Lind_max(1)-1;
    Lind_min = [Lind_min1 Lind_min2];
    lat_min = in.Wtime(Lind_min);
    amp_min = Wlp_3K(Lind_min);
    amp_min_bp = Wbp(Lind_min);

    %lat_max, lat_min, amp_min, amp_max_bp, amp_min_bp, amp_peak_470, lat_peak_470, lat_peak_470_last, amp_peak_470_last
    o.lat_max = lat_max;
    o.lat_min = lat_min;
    o.amp_min = amp_min;
    o.amp_max_bp = amp_max_bp;
    o.amp_min_bp = amp_min_bp;
    o.amp_peak_470 = amp_peak_470;
    o.lat_peak_470 = lat_peak_470;
    o.amp_peak_470_last = amp_peak_470_last;
    o.lat_peak_470_last = lat_peak_470_last;
    o.Wbp = Wbp;
    o.Wlp_3K = Wlp_3K;
    o.Wlp_470 = Wlp_470;

    % Calculate derived values from ECochG points for metrics export (o2)
    o2.SN = in.SN;
    o2.SPstar_amp = o.amp_peak_470(1,1); % SP* amplitude
    o2.Cstar_amp = o.amp_peak_470_last; % C* amplitude
    o2.SPstar_lat = o.lat_peak_470; % SP* latency
    o2.Cstar_lat = o.lat_peak_470_last; % C* latency
    o2.APstar_amp = o.amp_max_bp(1) - o.amp_min_bp(1); % AP* amplitude
    o2.N1P1star_amp = o.amp_max_bp(1) - o.amp_min_bp(2); % N1P1* amplitude
    o2.N1divN2 = o.amp_max_bp(1) / o.amp_max_bp(2);
    o2.N1divN3 = o.amp_max_bp(1) / o.amp_max_bp(3);
    o2.N1divN5 = o.amp_max_bp(1) / o.amp_max_bp(4);
    %o2.N1star_amp = o.amp_max_bp(1,1) - o.amp_min_bp(1,1); % N1* amplitude
    o2.N1star_amp = o.amp_max_bp(1);
    o2.N2star_amp = o.amp_max_bp(2);
    o2.N3star_amp = o.amp_max_bp(3);
    o2.N5star_amp = o.amp_max_bp(4);
    o2.N1star_lat = o.lat_max(1);
    o2.N2star_lat = o.lat_max(2);
    o2.N3star_lat = o.lat_max(3);
    o2.N5star_lat = o.lat_max(4);
    o2.N5starN1star_lat = o.lat_max(4) - o.lat_max(1);
    % close(hfig)
%DB
end